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Abstract. We analyze numerically ensembles of tight-binding Hamiltonians describ- 
ing highly-symmetric graphene nanoflakes with weak diagonal disorder induced by 
random electrostatic potential landscapes. When increasing the disorder strength, 
statistical distribution of energy levels evolves from Poissonian to Wigner, indicating 
the transition to quantum chaos. Power laws with the universal exponent map the 
disorder strength in nanoflakes of different sizes, boundaries, and microscopic disorder 
types onto a single parameter in additive random-matrix model. 
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1 Introduction 

Soon after the discovery of graphene — an atomically-thin monolayer of carbon 
atoms arranged in a honeycomb lattice llj — it was shown experimentally that 
electrons in this material behaves as spin-1/2 massless Dirac particles [2], in 
agreement with much earlier theoretical prediction by Semenoff For this 
reason, the nanostructures in graphene have attracted much attention, lead- 
ing physicists to reexamine classic effects of quantum transport [4 in search 
of novel features that arise from the unusual conical band structure, chirality, 
or the presence of additional quantum number [valley index) |5I6| . In par- 
ticular, a Coulomb-blockade experiment on quantum dots consist of graphene 
nanoflakes and normal metallic leads [7 shown signatures of quantum chaos 
(the energy-level repulsion) for the flake size smaller then 100 nm, but without 
clear identification of the system symmetry class. Some more light was shed 
on this issue with theoretical work [8 , showing that measurable quantities 
may indicate different symmetry class in the case of open than closed quantum 
dot. Later, the energy-level statistics of closed and irregular graphene flakes 
obtained from numerical diagonalization of tight-binding Hamiltonians |9ll0j 
was found to coincide with those given by the Gaussian orthogonal ensemble 
(GOE) of random matrices pTI] . 

In this paper, we follow the numerical approach established by Refs. |8l9)10j 
but focus on regular (hexagonal) graphene flakes with a weak diagonal disorder 
attributed to the substrate-induced random electrostatic potential landscape 
(see Fig. [I]). The results show, that the energy- level statistics of such systems 
coincide with those given by additive random matrices of the form H'^ + XV 
[12] . where i?" is the diagonal random matrix (and thus has Poisson statistics) 
and V is GOE matrix. We also found, that the parameter A is related to the 
extensive quantity A'totATg (where Atot is the total number of carbon atoms 




Fig. 1. Schematic representation of tlie systems studied numerically, (a), (b) Hexa- 
gonal graphene nanoflakes with armchair and zigzag edges with their radii Ra, Rz- 
(c) Conical dispersion relation E{kx,ky) near the Dirac point, (d) Typical potential 
profiles along the flake. Shaded areas on panels (c), (d) mark the energy range used 
when discussing the spectral statistics (see the main text for details). 

and Kq is an intensive measure of the disorder strength) via the scahng law 
A (X (A^tot^o)") with a ~ 0.6 regardless boundary conditions and microscopic 
details of the disorder model. 

The paper is organised as follows. In Sec. [2] we recall the basic findings on 
possible symmetry classes of chaotic nanosystems containing Dirac fermions, 
and present microscopic models of disorder in graphene nanoflakes. In Sec. [3j 
the random matrix model describing the transitions to quantum chaos is applied 
to rationalize level-spacing distributions obtained from numerical diagonaliza- 
tion of tight-binding Hamiltonians. The conclusions are given in Sec. |4] 

2 Dirac fermions in disordered graphene 

In this Section, we present two different microscopic models of disorder in 
graphene nanoflakes, representing the random electrostatic potential landscape 
abruptly or smoothly varying on the length-scale of the lattice spacing a = 
0.246 nm. But first, let us briefly recall (after Ref. [5]) the discussion of possible 
symmetry classes of such nanosystems. 

2.1 Symmetries of the Hamiltonian 

The effective Hamiltonian for low-energy excitations of electrons in graphene 
in the absence of magnetic field has a form of the Dirac Hamiltonian 

"Hoff = vfPxCTx 'S)Tz+ VFPyO-y «) To + [M{x,y)az + U{x,y)ao] (E) tq, (1) 

where vp w lO^m/s is the energy-independent Fermi velocity, <7i and {i = 
1,2,3) are the Pauli matrices acting on sublattice and valley degrees of freedom 



(respectively), and (Tq (''o) denotes the unit matrix. M{x,y) and U{x,y) are 
the mass term and the external electrostatic potential. Symmetries of the 
Hamiltonian ([T]) are defined by the following antiunitary operations: standard 
time reversal T, and two "special time reversals" 

T ^ {<JQ (E) Tx)C , Ts\^ ~i{(Jy<E)To)C, Tv ^ -i{(Jo<E)Ty)C, (2) 

where C denotes complex conjugation. The mass term breaks the symplectic 
symmetry associated with 7^i, leading to the two distinct possible scenarios: 

(i) In the case of weak intervalley scattering, Ty commutes with HcS, so 
the system consists of two independent subsystems (one for each valley) . Each 
subsystem lacks time- reversal symmetry, as T commutes only with full HcS- 
Because the Kramer's degeneracy {T^ — —I), the Hamiltonian consists of two 
degenerate blocks, each of which belonging to the Gaussian Unitary Ensemble 
(GUE). The analogous scenario was considered by Berry and Mondragon [13] 
for neutrino billiards, lacking the valley degrees of freedom. 

(ii) In the case of strong intervalley scattering caused by irregular and 
abrupt system edges, or by the potential abruptly varying on the scale of 
atomic separation, the two sublattices are also nonequivalent, so both spe- 
cial time-reversal symmetries T^i and Ty became irrelevant. T commutes with 
"HeS leading to the orthogonal symmetry class. 

The existing numerical studies for closed systems of irregular shapes |8I9I10| 
show that the typical intervalley scattering time is always shorter than the time 
required to resolve a level spacing (Heisenberg's time) leading to the scenario 
(ii). Some features of the scenario (i) were found in open systems fH^, for which 
the intervalley scattering time needs to be compared with much shorter escape 
time. Such systems are, however, beyond the scope of this paper, as we focus 
on regular and weakly-disordered systems, for which the intervalley scattering 
itself may be suppressed. 

2.2 Disorder in the tight-binding model of graphene 

The lattice Hamiltonian for disordered graphene reads 

n = E^'^-I*)01 + E [t^gatc(r.) + ;7i„,p(r,)] |z)(z|. (3) 

ij i 

The hopping-matrix element %j = —7 if the orbitals \i) and \j) are nearest 
neighbors on the honeycomb lattice (with 7 = ^^/Shvp/a. ~ 3eV), otherwise 
7ij = 0. The electrostatic potential contains a contribution [/gate from gate 
electrodes (slowly varying with the site position r^) and a random contribution 
Uimp from impurities. For small nanoflakes one can choose C/gato — ^^o = 0, 
whereas a realization of disorder potential is generated by randomly choosing 
-^imp lattice sites R„ {n ~ 1,. . . , Mmp) out of iVtot, and by randomly choosing 
the amplitudes C/„ S {—6,6). The potential is then smoothed over a distance ^ 
by convolution with a Gaussian, namely 



The special case of ^ <C a, A^imp = -^tot corresponds to the Anderson model on 
a honeycomb lattice, considered in work [T3] on spectral statistics of nanotube- 
like structures. Earlier, the model constituted by Eqs. (3j4| with ^ :$> a 



was shown to reproduce basic transport properties of disordered mesoscopic 
graphene samples |15|16) . It has not been considered, however, in the discus- 
sion of spectral statistics of nanoflakes so far. 

We further define the Fourier transform of two-point correlation function 

^ Ntat Wtot 

^1 " 777^^ — ^2 y] y] {Uhnp{n)Ui^piri)) exp [iq ■ (r; - r^)] , (5) 

(NtotflVF) 

where the system area A — |-\/3iVtota^, and the averaging takes place over 
possible realizations of the disorder Q (so (Ui^p) = 0). For the length scales 
large compared to ^, the dimensionless correlator 

^ y3iVi,,p f6Y 2 Jl, ife«a, ... 



9 A^tot \lj ' \|v^7r(e/a)2, if ^ » a, 

becomes a representative measure of the disorder strength. For q 7^ 0, we 
obtain Kq = Kq if a, or = i^o exp(— q^^^) if ^ >• a. The numerical 
value of the ratio Kq/Ko at q = (±|^,0) approximates the intervalley scat- 
tering rate, and is as small as 2 x 10~^ for ^ = VSa (used in the numerical 
simulations presented in remaining parts of the paper) . 



3 Random matrices and spectral statistics 

3.1 Additive matrix model for transition Poisson-GOE 



Before presenting the numerical results for spectral statistics of graphene nanoflakes, 
let us briefly review corresponding additive random-matrix models and result- 
ing nearest-neighbor spacings distributions ^12) . 

When large integrable system undergoes transition to quantum chaos, its 
spectral properties can be modelled by the following random Hamiltonian 



H = 



H° + XV 



(7) 



where i?" is diagonal random matrix, which elements follow a Gaussian dis- 
tribution with zero mean and the variance {{Hfj)'^) = Sij, the parameter 
As [ 0, 00 ] , and 1^ is a member of one of the Gaussian ensembles. In particular, 
for transition Poisson-GOE, elements of V are real numbers chosen to follow 
a Gaussian distribution with zero mean and the variance (Vj^) = (1 -I- Sij)/N, 
where A^ is the matrix size. 

For N — 2, the eigenvalue-spacings distribution for the Hamiltonian Q can 
be found analytically and reads, for transition Poisson-GOE, 
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Fig. 2. Nearest-neighbor spacing distribution P{S) for hexagonal flakes of Fig. [l] 
(a)-(c) Armchair edges, Anderson model of disorder = 0, Nimp ~ Ntot = 8322). 



(d)-(f) Zigzag edges, smooth impurity potential = \/3i, Mmp ^ Nu 



10584). 



Disorder strength Ko ([gJi is varied between the panels by changing the potential ampli- 
tude 6 [panels (a)-(c)] or by fixing 5/^ — 0.1 and varying the impurity concentration 
A'imp/A'tot [panels (d)-(f)]. Histograms show the numerical data obtained by aver- 
aging over 200-400 disorder realizations. Solid lines show the statistics interpolating 
between Poisson and GOE ([sjl with best-fitted parameter A — Aflt specified for each 
panel. The limiting cases of Poisson (A = 0) and GOE (A = oo) statistics are shown 
with dashed and dotted lines (respectively). 



Io{x) is the modified Bessel function of the first kind; u{X) = y7rf/(— |, 0, A^) 
with U{a,b,x) the confluent hypergeometric function |17j . In particular, for 
A = the Poissonian distribution P{S) = exp(— 5') is restored. For the opposite 
hmit (A — > oo) we have P{S) — (7r/2)S'exp(— 7rS'^/4), reproducing the Wigner 
surmise for GOE matrices. For < A < oo, Eq. Q describe level-spacings 
distributions interpolating between Poisson and GOE statistics, with P(A; S) oc 
S^/A if 5 < A < 1, or P(A; S") cx 5 if < 1 < A. 

For large N, the statistics P{S) (so-called nearest-neighbor spacings distri- 
bution) is defined as a distribution of a variable S = (£"«+! — En){p), where 
(p) is the average density of states, and £"„ < En+i are neighboring energy 
levels. Subsequently, we have P{S) = SP{S) — 1 (so-called unfolded 
spectrum) . Although Eq. ([8| is exact for = 2 only, it was shown numerically 
[12j that P(Afit;S') with Afit — \/NX provides an excellent approximation of 
P{S) for large random matrices of the form given by Eq. (It]). 




Fig. 3. Least-squares fitted parameters Aflt for transition Poisson-GOE ([Sf as func- 
tions of disorder strength for hexagons with armchair edges (a) [or zigzag edges (b)], 
different sizes, and the two distinct disorder types. A'tot ~ 8322 (a) [or 10584 (b)]: 
f = (□), C = (O); Ntot = 34062 (a) [or 42366 (b)]: ^ = ^/3a (v); 

Ntot = 6144 [panel (b) only]: ^ = (■), ^ = VSa (a). Lines denote best fitted 
power-law relations for the two disorder types (see Table [T] for details). 



3.2 Energy- level distributions for disordered graphene flakes 



In this Subsection, the central question of the present work is addressed, 
namely: Whether the statistic interpolating between Poisson and GOE, P(A; S) 
(|8| is capable of describing nearest-neighbor spacings distributions P{S) for 
weakly-disordered graphene flakes? In other words, may the additive-random 
matrix model defined via Eq. ^ he applicable for such relativistic nanosys- 
tems? To answer this question, we focus on two systems of a high symmetry: 
hexagonal flakes with entirely armchair or zigzag edges, each of which is show- 
ing Poisson statistic in the absence of disorder (providing the level degeneracy 
is properly taken into account). As already mentioned in Sec. [2j two distinct 
models of disorder are applied to eac h sy stem: Anderson model, defined by 
setting ^ = and Mmp = A'tot in Eqs. (3 4), or smooth disorder, with ^ = \/3a 
and Aimp < A^tot- 

To obtain the statistics P{S), we diagonalized numerically tight- binding 
Hamiltonians ^ for the flake containing A^tot S 10'' atoms and 200 — 400 
independent disorder realizations for either type of edges, disorder models, and 
each disorder strength quantified by the correlator Kq ([6|. Some additional 
effort is required when unfolding the spectra: Unlike for two-dimensional gas 
of Schrodinger electrons, for which average density of states (p) is assumed to be 
energy-independent, for bulk graphene we have [TSj {p{E)) ~ ^|£'|/[7r(/iw^)^] 
(per spin). For small systems studied here, boundary effects lead to additional 
states appearing near E ~ ±7 (armchair edges) or i? ~ (zigzag edges). 
Also, the impurity potential ^ introduces some bound states for \E\ < S. All 
these additional states, however, are localized on areas small in comparison to 
A, and thus not contribute to the spectrum obtained in a Coulomb-blockade 
experiment such as reported in Ref. 7J. For this reason, we limit the energy 



Table 1. Least-square fitted power-laws Afit(C) ~ -^iC"j with ^ = Nt^tKo (lines in 
Fig. [3|. Numbers in parenthesis are standard deviations for the last digit. 



Disorder model 


Armchair edges 


Zigzag 


edges 


C=0, iVimp^iVtot 
^=v/3o, Afimp<Artot 


Ai= 0.059(2) a= 0.55(1) 
0.035(2) 0.56(1) 


Ai= 0.046(3) 
0.023(4) 


a= 0.59(1) 
0.56(3) 



range [cf. i?min and -Emax in Fig. [TJc),(d)] such that 

{p{E))c,po + lj^\E\, for 0.1 < |i?|/7 < 0.5. (9) 

The constant term po and the effective area AcS ^ A. are determined via least- 
square fitting of Eq. ([9| to the actual {p{E)) obtained by numerical averaging 
over independent disorder realizations. 

Our numerical results are presented in Figs. [2] and ^ First, we compare 
the statistics P{S) on two selected examples of nanosystems considered: the 
hexagon with armchair edges and Anderson- type disorder (Fig. [2]ja)-(c)) and 
the hexagon with zigzag edges and smooth disorder (Fig. [2jd)~(f)). Although 
some systematic deviations of P{S) from the best-fitted interpolating statistics 
P(Afit;5') (|8| are visible for S > 1 due to a finite system size (notice that a 
better agreement is observed for iVtot = 10584 than for 8322), P[\f^t] S) repro- 
duces the actual nearest-neighbor spacings distribution with a good accuracy 
for both systems and wide range of Kq. We further notice, that similar values 
of Afit are reached for the second system at K^) typically 5 — 8 times larger than 
for the first system. 

The dependence of Afit on the total disorder strength NtotKa for all datasets 
available is illustrated in Fig. [s] (datapoints) in the logarithmic scale. The 
particular choice of the independent variable C = A^tot^o allows us to find the 
approximating relations A^t ~ Afit(C)7 which still differ between the systems 
with different edges or disorder types, but remain unchanged when varying A^tot 
and Kq independently with the remaining parameters fixed. Also, for smooth 
disorder, we vary A'imp having 5 fixed at 6/^ = 0.1 or 0.5 (corresponding to 
the absence or presence of charge puddles in the physical system) . 

Least-square fitted power-laws Afit(C) are listed in Table [l] and plot in Fig. 
[3][a),(b) (lines). The power-laws fail for A^t ^ 1, as P(Afit; S) becomes indistin- 
guishable from GOE statistics in this range. They are, however, closely-followed 
by the datapoints for smaller Agt-s. We further verify the obtained Afit(C)-s for 
the case of smooth impurity potential, by taking the systems approximately 
four times larger in area (namely, A'tot = 34062 for armchair edges or 42366 
for zigzag edges), but generating only one disorder realization for each Kq. 
Such an approach reproduces the experimental procedure of Ref. [7:, where 
the spectrum of a single system was obtained. Additionally, the corresponding 
flake diameters 2Ra = 87^/3 a ~ 37 nm and 2Rz = 168 a ~ 41 nm are of the 
same order of magnitude as diameters reported in Ref. [7^ . The new datapoints 
(open triangles in Fig. [3]) still follow the corresponding power-laws, providing 
that Afit < 1. 



Probably, the most remarkable feature of these results is that all graphene 
nanoflakes considered show transition Poisson-GOE when increasing the dis- 
order strength, with no signatures of GUE statistics. This is expected for the 
flakes with armchair edges which couple the valleys [S], or with zigzag edges 
and Anderson- type disorder [15], for which the intervalley scattering restores 
time-reversal symmetry. The absence of GUE statistics seems surprising in 
the case of zigzag edges accompanied by the smooth impurity potential. In 
such case, some intervalley scattering originates from six 120° corners, a role 
of which may become decisive for spectral statistics of closed nanosystems. 

4 Conclusions 

We find that the additive random-matrix model, describing a transition to 
quantum chaos in Hamiltonian systems, is also relevant when discussing spec- 
tral statistics of highly-symmetric graphene nanoflakes with a weak diagonal 
disorder. The functional relation between the model parameter A and the dis- 
order strength N^otRQ has a form of a power law, with the universal exponent 
a ~ 0.6, which is insensitive to the boundary type or to the microscopic model 
of the impurity potential. 

In the chaotic range, regular graphene flakes show energy-level statistics 
characteristic for the Gaussian Orthogonal Ensemble (GOE) of random ma- 
trices, indicating the strong scattering of Dirac fermions between the valleys. 
This coincides with earlier findings for irregular nanoflakes |8I9I10) . 
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